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We present a new method for computing the transverse transfer matrix through superimposed ax- 
isymmetric RF and solenoid field maps. The algorithm constructs the transfer matrix directly from 
one dimensional RF and solenoid field maps without computing numerical derivatives or eigenfunc- 
tion expansions of the field map data. In addition, this method accurately describes the dynamics 
of low energy particles starting from a solenoid immersed cathode, allowing the method to be used 
to simulate transport through both RF and electrostatic guns. Comparison of particle tracking with 
the transfer matrix and direct integration of the equations of motion through several field set-ups 
shows excellent agreement between the two methods. 



I. INTRODUCTION 

Linear transfer matrices continue to serve as important simulation tools in the design, commissioning, and operation 
of modern accelerators. Common examples of their use include computing linear centroid motion, beam based 
alignment of optical elements, and orbit feedback. In addition, they also feature prominently in the theory of round- 
to- flat beam transforms in RF and DC guns .1, 2], and are used for emittance measurements for beams without space 
charge [3]. Because of this utility, analytic expressions for the transfer matrix through many beam line elements, such 
as magnets with constant fields, are well known and are widely used. In contrast to these simple elements, the fields of 
many beam line elements in modern accelerators have no analytic form and may overlap each other. For example, in 
high brightness electron sources, solenoid fields used for emittance compensation may overlap the accelerating fields at 
or near the cathode. To properly describe the dynamics in these machines, the transfer matrix through superimposed 
RF and solenoid fields must be constructed. In general, to model these elements one must use numerically computed 
electromagnetic field maps. Unfortunately, no closed form solutions for the transfer matrix through such elements 
exist. 

Nonetheless, a significant amount of work has gone into developing both semi-analytic and numerical techniques 
to compute these matrices. In general, these techniques require some form of manipulation of the field map data 
and often have a limited range of validity. For example, the widely used RF transfer matrix given by Rozensweig 
and Serafini [4] requires a Floquet expansion of the on-axis RF field map, and is only valid for ultra-relativistic 
particles. Other methods expand the field map data in terms of the the general solution to the homogeneous Maxwell 
equations in cylindrical coordinates [SlIH]- From this expansion, the vector potential can be computed, allowing the 
Hamiltonian for the overlapping solenoid and RF fields to be constructed. When combined with differential and lie 
algebra techniques this method is quite powerful, and can be used to generate arbitrary order maps. Despite this, it 
requires significant overhead in setting up and may not be suitable for online modeling purposes. A simpler solution 
has been put forth in [7] . In this approach, the particle energy and RF fields are assumed constant over a small time 
step. The second order linear transverse equations of motion can then be solved 'exactly', and the transfer matrix for 
the time step constructed. While this method works well for relativistic particles, it does have several draw backs. 
First, the assumption that the energy is constant requires extremely small time steps for very low energy particles 
like those emitted from a photocathode. Secondly, the determinant of the resulting transfer matrix is only correct to 
first order in the step size. 

Depending on the particular application, one or more of these methods may be appropriate to use. However, 
the need for a simple algorithm that computes the transfer matrix through superimposed RF and solenoids fields 
accurately for low particle energies still exists. In order to be useful, such a matrix should satisfy the following three 
requirements: (i) using reasonable step sizes, the matrix should be able to describe the low energy dynamics found in 
RF or DC guns, (ii) the matrix should have the correct determinant regardless of the step size, (iii) the algorithm for 
constructing the matrix should be easy to implement. Based on a new solution to the equations of motion, we derive 
a transfer matrix with all three of these qualities. 

The layout of this work is as follows. First, the longitudinal equations of motion are solved for a single small step in 
the independent variable. Then, the transfer matrix over the same step is derived for electrostatic and solenoid fields. 
Building on this result, the matrix for combined RF and solenoid fields is computed. This matrix is then tested with 
tracking through a DC gun, a superconducting RF cavity and a RF gun with a solenoid immersed cathode. Excellent 
agreement between the transfer matrix and direct integration of the equations of motion is demonstrated in all three 
cases. 



II. FIELD EXPANSIONS AND THE EQUATIONS OF MOTION 

In general, both standing and traveling wave RF fields can be written in the complex form E = £(a;, j/,z)e'"* and 
B = B{x, y, 2)6*"^*. In these and all subsequent expression, tildes are used to denote phasor quantities. The coordinate 
system for the fields is set up so that the z-axis points along the length of the beamline. The functions £ and B 
represent the field maps generated by RF cavity field solvers, and include the initial phase offset of the cavity. For 
notational simplicity, the phase factor e*"* and the real symbol Re[...] are suppressed in all but the final results for 
the transfer matrix. To derive the linear equations of motion, the fields are expanded to first order in the transverse 
offsets X and y. The linearized RF fields take the form 




^>.y--) --7^V^]^-U^]y + £^ir = o)z 



y n-'^\ d - , ^ n^ 



m->y,-) = -|(9^^x+9P:^)£.y. (1) 



Similarly, the expansion of the solenoid field gives 



„ X / dBz \ v ( dBz \ „ , ^ 

Note that the solenoid field is distinguished from the RF magnetic field by the lack of a tilde. The total physical fields 
are given by Etot = Re[f e*"*] and Btot = Re[^e'"*] + B.oi- 

For the single particle equations of motion, the longitudinal coordinate z is used as the independent variable. 
Derivatives with respect to z are denoted with a prime: /' = df /dz. The equations of motion for the longitudinal 
phase space variables t and 7 are 

t' = ^, f(.) = ||=^V^e, y = Re[7'e-]. (3) 

Here the constant £e = mc^ /e gives the (signed) rest energy of the electron in [eV]. The normalized energy is given 
by 7 = Re [76*"*]. The transverse equations of motion can be written as [Hj: 

2AeW~(^)y^o, (4) 




28eP 

2A9'^x'+(^]x = 0. (5) 



2£eP 

In this expression, p ~ /3j is the normalized reference particle momentum, and A9l is the Larmor angle defined by 

AeL{z) = - r ^dz, Ae'L = -^. (6) 

Note the use of a negative sign in front of the integral. With this definition, a positive solenoid field creates a positive 
change in the Larmor angle for electrons. Decoupling these equations requires rotating to the Larmor frame. By 
defining the variable rj = x -\- iy, the transformation to Larmor coordinates takes the form [8] : 

VL{z)^XL+iyL^Ve-''^'-^'l (7) 

The equivalent transformation for the transverse phase space vector u — [ x x' y y' j is defined by 
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In this and following expressions, C = cosA^/,, and S = sinA^^. The matrix in this equation is denoted by 
L — L{AO L , AO'i^) . It is easy to show that the determinant of this matrix is equal to unity. Another important 
characteristic of this matrix is found by taking the limit as z — > z^. In this limit, AO^ — > 0, and the matrix reduces to 
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The two remaining focusing terms describe the effect of starting a particle with B^ 7^ 0, the so called 'immersed 
cathode' condition. 

Substituting the transformation in Eq. (l7| into Eqs. Q and (Is]) gives the transverse equation of motion for the 
Larmor coordinates: 



p'^ \ 2£,p2 ^ 2cp 



^U ^.U Ufz. + ^ + (A^l)M .L = 0. (10) 



In this expression r]L stands for either xl and yi,. Because the transformation between the laboratory and Larmor 
coordinates is known, all that remains is to solve the above differential equation. With this solution, the transfer 
matrix for the lab phase space variables follows directly: M{zi — > zj) = L~^{zf)Mij{zi -> Zf)L{zi). Consequently, 



the majority of this work is spent finding an exact solution to Eq. ( 10 1 within the approximation that the fields are 
constant over a small step in z. 



III. DERIVATION OF THE TRANSFER MATRIX 



As stated before, general solutions to the differential equations in Eqs. (SlIS), and (10 1 do not exist for arbitrary 



cavity and solenoid field maps. As a result, any method for computing the transfer matrix analytically requires 
some form of approximation to these equations. The approach taken in this work is to find an exact solution to the 
equations of motion for a step in z and t that is small enough so that the fields profiles E-ziz) and Bziz), as well as 
the RF phase, don't change appreciably. The solution is exact in the sense that the particle energy changes correctly 
over the course of the step. The change in the field map profile as well as the RF phase are then included with the 
use of edge focusing matrices. Slicing the field maps and consecutively multiplying the matrices for each step gives 
the total transfer matrix: 

M(z, ^z/)« J|AM(zfc^Zfc + Az). (11) 

fe 

The first step in constructing the transfer matrix for one step is to solve the longitudinal equations of motion in 
Eq. ([3]). For a constant electric field, 7' is constant, and the normalized energy, momentum, and velocity are given by 

7(z)=7,+7'(z-z,), p(z) = V(7»+7'(^-2.))'-l, /3(z)=p(z)/7(z). (12) 

Using these expressions, the derivatives of i(z) and A6'l(z), in Eqs. (l3| and (10) can be directly integrated: 



t(z) = i-b(z)-p.], (13) 

A..W^(A),„(?«±l«)^ (14) 

In the last line, the constant h is defined as h = p ■ A0^ = —eBz/2mc. In addition to defining the transformation 
between the lab and Larmor coordinates, the function A6'i plays an important role in the derivation of the transfer 
matrix for both the electrostatic and RF field cases. 

A. Overlapping Electrostatic and Solenoid Fields 

The derivation of the transfer matrix for superimposed RF and solenoid fields is based in part on the method used 
to derive the transfer matrix for static fields ^ . In this section, a detailed derivation of the static field result is given. 
The same techniques are then modified and used to derive the RF matrix in the following section. For static fields. 



the equation of motion is found by taking w ^- in Eq. ( 10 ) 



Note that p' ex Ez and A0^ oc Bz in this expression, implying "q^ depends on both the accelerating field and its 
gradient, as well as the solenoid field. The transfer matrix from Zi to Zf = Zi + Az is derived in a three step process. 
Over the interval [zi^Zf], the electric and magnetic fields are approximated as rectangular step functions. Formally 
the fields are written as: 
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B', = Bz{zi){5{z-Zi)-5{z~Zf)}, 



(16) 

where 9{z) is the Heaviside step function and 5{z) is the Dirac delta function. The transfer matrix is then found by 
solving the transverse equation of motion piecewise from Zi and Zf . 

First, the equations of motion are integrated across the rising edge of the electric field at Zi. Because the rising 
edge is approximated as an instantaneous step, the particle's position does not changed: rj^^zf) = riL{z~) = rj^izi). 
Integrating the equation of motion gives the kick delivered to the particle's trajectory: 



^i'l = - t-?;, + U& + (A«i)M '« M" - / i^T^'^i-i'W' -'■)''= = -7d^iU'i) 



The corresponding transfer matrix for the rising edge, Re, takes the form 

REh,l')=( _\_ °V (17) 

\ 27/32 -L / 



Next, the equation of motion is solved across the interval {zi,Zf) where both the electric and magnetic field are 
approximately constant. In this region, the equation of motion reduces to 
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(18) 



where A6'^(z) ex l/p{z). With the electric field held constant, 7, p, /3, and A9l are given by Eqs. (12) and (14). With 
these functions, the differential equation can be solved by assuming rj^ = 77L(A6'i). Plugging this into Eq. (18 1 gives 
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(19) 



Using Eq. (14 1, is is possible to show IS.9" — — 7'A6''/7/3 , canceling the second term the above equation. Assuming 
AS' 7^ 0, the first term in this expression must also vanish. It follows that: 77^ = Acos/S.Ol + -BsinA6'i. Completing 
the initial value problem for this solution determines the transfer matrix for the step from Zi to Zf: 
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(20) 



The last step in constructing the full matrix for the interval [zi , Zf] is to evaluate the transfer matrix for the falling 
edge of the accelerating field. The result is essentially the same as before, except now the derivative of the electric 
field has the opposite sign. This allows the transfer matrix for the falling edge to be written as 
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7' 

27/32 



Combining the three matrices for each region gives the full transfer matrix for the step Az: 



(21) 
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(22) 



One important thing to note about the matrix in Eq. ( 22 ) is that it has the correct determinant for the phase space 
variables chosen: det(AM^';j.,) = Pi/pf- The transfer matrix for the canonical phase space variables xl and Px,L can 
be found by applying the transformation: 
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dc 







(23) 



It follows from this expression that the matrix AAf^p^ satisfies the symplectic condition det(AM^p^) = 1. In addition 
to this, the transfer matrix also has the convenient feature that the derivative of the accelerating field never has to 
be calculated, bypassing the need to compute derivatives numerically. It is important to note that when starting a 
particle from a non-zero electric field, the first rising edge matrix should not be included . Doing so amounts to the 
particle seeing the field rise from zero to the actual value of the field at the cathode and is not physical. Similarly, 
the falling edge matrix at the end of tracking should not be included when the particle is found in a non-zero electric 
field. 



B. Overlapping RF and Solenoid Fields 



With the results for the electrostatic field worked out, it is now possible to construct the transfer matrix for RF 
fields. The approximation used here is similar to that used in the electrostatic case: the field profiles £z and B^, 
as well as the RF phase, are assumed constant over the step Az. This implies j{z,t) is also constant, allowing the 
solutions to the longitudinal equations of motion in Eq. ( 12 ) to be used. The effect of adding the time dependence to 



the edge matrices is minimal; the rising edge matrix remains the same. For the falling edge matrix, the change in the 
RF phase is included in the electric field: j'f = 7'e*"'^*. Note that this reduces to the electrostatic case when a; ^- 0. 
The equation of motion over the interval (z^, zj) reduces to 



^^ + 1^^^ 
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(24) 



Unfortunately, the extra focusing from the RF magnetic field scales as p~^, and the form of the solution used in the 
electrostatic case no longer solves the equation of motion. To our knowledge, there is no solution to this equation in 
its current form. 

In order to solve this equation, it must be transformed in such a way that the RF magnetic focusing term is removed. 
To do so requires switching the independent variable from longitudinal position to time. The resulting transformation 
of the phase space variables takes the form 
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The transverse equation of motion with time as the independent variable is given by 
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The second transformation is to reduce this equation by introducing the variable ^l = s/^tjl. 
transformation can be written in matrix form: 



7L 



= A 






/ 1 

\ 27 



1 



The equation of motion for the reduced variables takes the form of Hill's equation: 



i) + 



(p2 



4(p2 



1 



+ 



2j^£o 



{Ml 



fj^O. 



(25) 

(26) 
Once again, this 

(27) 
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The important thing to notice about this expression is that the problematic RF focusing term has been hidden in 
the variable transformation. All that remains now is to solve this equation in the region where the solenoid and 
accelerating fields a re c onstant. To do so requires knowing the functions p{t) and 7(i). The momentum is easily found 
by rearranging Eq. ( 13 ): p(t) = pi + c~f' {t — ti). The normalized energy is then given by 7(i) = ^Jp'^{t) + 1. Inserting 



these expressions into the equations of motion in the constant field region yields: 
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This equation can be solved with the function ul = yf^{A cos 1^6 l + Bsin/S.OL). The transfer matrix is found by 
completing the initial value problem for this solution. The resulting matrix can be written in the compact form 



Mt,f = A(7/,^/)T(/3;)A/d4^r-i(ft)A-^(7.,7). 



(30) 



In this expression 7f = c/3f^' . This matrix correctly describes the evolution of the reduced variables. In order to get 
the transfer matrix for the usual phase space variables it must be transformed back: 



M^l^j = T-'i(3f)A-\jf,^^e^-^')A{^f,%)Til3f)Mf^ 



(31) 



It is clear from this expression that the effects of the RF magnetic focusing must be contained in the matrices left 
multiplying Mf_^r. Combining these matrices together gives 
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In the limit that At = t — tj is small, this matrix reduces to 
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From equation Eq. (10 1, it is clear that this is nothing but a thin lens approximation to the focusing delivered by the 
RF magnetic field. The complete matrix for the step Az is found by including the effects of the field edges: 

AM:%, = R-^' (7/,7'e*"^*) i?rf(7/,7 ,e''^^*)Mf4/i?i.(7.,7 ) 
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(34) 



This is the main result of this work. As a reminder, 7' = e£ z{zi) / rnc? , where £z is the complex electric field map. 
Additionally, the expressions for 7, p, and At can be found in Eqs. ( 12p4 ), respectively. For clarity we leave the 
result in the above factorized matrix form. This allows several limiting cases to be easily evaluated. First, in the 
limit that both the RF and solenoid fields vanish, 6,7' -^ 0, and Eq. (34) reduces to a drift matrix. For vanishing RF 
fields, 7' — > 0, and Eq. (34) reduces to a hard edge solenoid matrix. In the limit that oj — > 0, i?rt reduces to /2x2, and 
the total RF matrix reduces to the previous electrostatic result in Eq. ( 22 ) . In addition to having the correct limiting 
behavior, the RF transfer matrix also has the correct determinant: det|AM^^^,] = det[AAf^'^^,] = VilVf- This follows 
directly from the fact that dct[i?if] = 1. 



IV. TESTING THE TRANSFER MATRIX 

To test the validity of our approach, the energy gain and transfer matrix are calculated through three different field 
set-ups and compared to direct integration of the equations of motion using a fourth order Runge-Kutta algorithm. 
The three field set-ups used are: a DC gun and with overlapping solenoid, a SRF cavity, and a RF gun with a solenoid 
immersed cathode. To check that the transfer matrix correctly describes the transverse dynamics in each case, all 
four transfer matrix elements are compared with Rugne-Kutta integration. To do so, the two principle trajectories 
through each field set-up are computed. These trajectories are defined by the initial phase space coordinates Ui = 
( 1 mm, ) , and U2 = ( 0, 1 mrad ) . Tablejl shows the various settings used in each simulation. 







TABLE I 


. Simulation 


Parameters 






Field Set-up 


Voltage 


Phase 


Max. Bsoi 


KE(z = 0) 


Step Size Type 


Step Size 


DC gun & Solenoid 


500 kV 





0.04 T 


1 eV 


fixed 


0.1 mm 


SRF Cavity 


3MV 


on-crest 


N/A 


1 MeV 


fixed 


2 mm 


RF gun & Solenoid 


1 MV 


on-crest 


0.04 T 


1 eV 


adaptable 


0.1 mm (avg.) 



For the first simulation, we use the field maps for the the high voltage DC gun and first omittance compensation 
solenoid of the Cornell ERL injector prototype. Fig. 1(a) shows the field maps corresponding to gun voltage and 



__ Fig. 
solenoid field strength given in Table-ll] Fig 1(b) shows 
in Eq 

field solution is 0.1 mm. The agreement between the two methods demonstrates that the constant field solution works 
very well for the longitudinal variables, even for low initial kinetic energies. Fig. [5] shows the results of tracking the 



low the energy gain computed from the constant field solution 
12 ) compares to the energy gain computed using Runge-Kutta integration. The step sized used for the constant 
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(a) Field maps for the dc gun and solenoid. 
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(b) Energy gain through the dc gun. 



FIG. 1. The fields and energy gain for a 500 kV gun voltage, 0.04 T maximum solenoid field setting, and 1 eV initial kinetic 
energy. 
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FIG. 2. Comparison of Runge-Kutta integration (blue) and the tracking using the transfer matrix (green) through the DC gun 
and solenoid fields. 



principle trajectories using both the RF transfer matrix with uj — 0, and Runge-Kutta integration. As with with 
longitudinal variables, the agreement between both methods of tracking is excellent. In addition to these results, 
the expression for the electrostatic transfer matrix has been experimentally verified in [3]. Next, the two principle 
trajectories are computed through the field map of the 1.3 GHz Cornell ERL injector SRF cavity. Fig. 3(a) shows 



the on-axis electric field map of the SRF cavity model with a 3 MV cavity voltage. The corresponding energy gain 
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(a) SRF cavity field map. 
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(b) Energy gain through the SRF cavity. 



FIG. 3. The field map for the SRF injector cavity and the corresponding energy gain. The cavity voltage is 3 MV, the initial 
kinetic energy is 1 MeV, and the phase is on-crest. 



through the cavity computed using the constant field solution and Runge-Kutta integration are shown in Fig. |3(b)[ 
The results of tracking the two principle trajectories through the cavity with a fixed step size of 2 mm are shown 
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FIG. 4. Comparison of direct integration (blue) and tracking using the transfer matrix (green) of the two principle trajectories. 



in Fig HI As in the electrostatic case, the agreement is very good. Finally, the principle trajectories are computed 
through the RF gun set-up. To simulate an RF gun, the last 1.5 cells of the injector cavity field map are used. The 
solenoid field is positioned so that the maximum value of the solenoid field occurs at the cathode. In order to make 
sure that the RF phase is constant over each step, a simple adaptive step size algorithm is included. This algorithm 
adjusts the step size so that the change in RF phase over the step is less than a user defined tolerance. Fig. 5(a) 
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(a) RF gun field maps. 
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(b) Energy gain through the RF gun. 



FIG. 5. The fields and energy gain . The cavity field is scaled and rotated so that the cavity voltage is 1 MV, and the phase 
is set to the on-crest value for a 1 MeV electron. 



shows field maps for the RF gun set-up. The corresponding energy gain through the gun is shown in Fig. |5(b)[ The 
accelerating field is scaled so that the RF cavity voltage is 1 MV and the phase is set for maximum acceleration. The 
tracking results for the two principle trajectories are shown in Fig. [6] From the figure, it is clear that that the transfer 
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FIG. 6. Comparison of direct integration (blue) and tracking using the transfer matrix (green) of the two principle trajectories, 
matrix works well in the low energy case. The average step size for the simulation was roughly Az = 0.1 mm. 

V. CONCLUSION 

We have derived and tested a new method for calculating the 4x4 transfer matrix through superimposed RF 
and solenoid fields. The algorithm computes the transfer matrix directly from the field data without computing 
eigenfunction expansions or numerical derivatives. Comparison with numerical integration demonstrates that this 
new method works for low energy beams starting from a solenoid immersed cathode. Additionally, because the 
algorithm relies on an analytic solutions to the equations of motion, it is simple to implement and guarantees the 
correct value for the determinant of the transfer matrix. One limitation to this approach is the assumption (inherent 
in the derivation) that the fields display cylindrical symmetry. For many applications this is a reasonable assumption, 
however previous work shows that asymmetric focusing from input power couplers may be noticeable when heavy 
beam loading is present [^. In addition, when tracking ultra-relativistic particles, the algorithm takes steps typically 
on the order of a few millimeters, and therefore may not be the best choice for computational speed. In this case, 
one may still choose to use the Rosenzweig-Serafini matrx. Nonetheless, the matrix algorithm given here strikes an 
appropriate balance between accuracy, speed, and simplicity not previously achieved. 
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